################################################################################
## Group Identities and Parliamentary Debates: Replication package
## Fiva, Nedregård and Øien (2025)

# Description:

## Code to make Figure A4: "Comparison of actual estimates and mean placebo estimates 
## for each election year"

################################################################################



# Packages

library(data.table)


# Directories (the parent directory is set by master.R)
data.dir <-  "../data/3_model_output"
fig.dir  <-  "../results/figures"

####################################

#---------GATHERING DATA----------#

####################################


data.list <- vector("list", length = 4)
names(data.list) <- c("gender", "age", "town", "occ")


charac <- c("gender", "age", "town", "occ")
labels <- c("Gender", "Age", "Urbanicity", "Social background")

for (c in charac){
  
  d   <- fread(file   = glue::glue("{data.dir}/{c}_comparty.csv"))
  
  d_p <- fread(file = glue::glue("{data.dir}/{c}_comparty_placebo.csv"))
  
  # Main estimates
  
  d[, type := "Main estimates"]
  
  d_p[, pi := rowMeans(.SD), .SDcols = names(d_p)[names(d_p) != "sessions"]]
  d_p <- d_p[, .(type = "Mean placebo estimates", pi)]
  d_p$session <- d$session
  
  d <- rbindlist(list(d, d_p), use.names = T)
  d[, mean_pi := mean(pi), by = "type"]
  
  d[, label := labels[which(charac %in% c)]]
  
  data.list[[c]] <- d
  
}

results <- rbindlist(data.list)
results[, label := factor(label, levels = c("Gender", "Age", "Urbanicity", 
                                            "Social background"))]

mean_results <- results[, .(mean_pi = mean(mean_pi)), by = .(label, type)]

####################################

#---------PLOTTING-----------------#

####################################

pdf(paste(fig.dir, "figA4.pdf", sep = "/"), width = 18, height = 15)

layout(matrix(c(1, 2, 3, 4, 5, 5), nrow = 3, byrow = TRUE), heights = c(1, 1, 0.2))

# Set up 2x2 plot layout
par(mar = c(5, 8, 2, 1), las = 1)

# Define the unique labels for the facets
labels <- unique(results$label)

for (lbl in labels) {
  # Subset data for each facet
  subset_data <- subset(results, label == lbl)
  
  # Plot the first histogram (Main estimates)
  hist.sum <- hist(subset_data$pi[subset_data$type == "Main estimates"], 
                   breaks = 40, plot = FALSE)
  hist.sum$density <- with(hist.sum, density * diff(breaks)[1])
  
  plot(hist.sum, freq = FALSE, xlab = "Estimates", xlim = c(0.5, 0.506), 
       ylim = c(0, .12), main = lbl, cex.main = 2.5, ylab = NULL,
       col = rgb(1, 1, 1, 0.5), cex.axis = 2, cex.lab = 2, border = FALSE)
  
  # Add the second histogram (Mean placebo estimates)
  hist.sum.placebo <- hist(subset_data$pi[subset_data$type == "Mean placebo estimates"], 
                           breaks = 40, plot = FALSE)
  hist.sum.placebo$density <- with(hist.sum.placebo, density * diff(breaks)[1])
  
  plot(hist.sum.placebo, freq = FALSE, xlab = "Estimates", xlim = c(0.5, 0.506), 
       ylim = c(0, .12),
       col = rgb(1, 1, 1, 0.5), add = TRUE, border = FALSE,
       cex.axis = 2, cex.lab = 2)
  
  # Add density lines for both types

  
  # Calculate density for Main estimates
  dens1 <- density(subset_data$pi[subset_data$type == "Main estimates"], adjust = 0.6)
  
  # Scale dens1 to match histogram density
  scaling_factor <- diff(hist.sum$breaks)[1]
  dens1$y <- dens1$y * scaling_factor
  
  # Add scaled density line for Main estimates
  lines(dens1, col = "gray40", lwd = 2)
  
  
  # Calculate density for Main estimates
  dens2 <- density(subset_data$pi[subset_data$type == "Mean placebo estimates"], adjust = 0.6)
  
  # Scale dens2 to match histogram density

  dens2$y <- dens2$y * scaling_factor
  
  # Add scaled density line for Main estimates
  lines(dens2, col = "black", lwd = 2, lty = 2)
  
  mtext("Density", side = 2, line = 5, cex = 2, las = 3)  
  
  
}



# Set up margins and parameters for the legend
par(mar = c(0, 0, 0, 0))  # No margins for the legend area
plot.new()  # Create a new, empty plot for the legend

# Add the legend
legend("center", 
       legend = c("Placebo", "Actual Estimates"), 
       lty = c(2, 1), 
       lwd = 2, 
       col = c("black", "gray40"), 
       horiz = TRUE, 
       bty = "n", 
       cex = 3)


dev.off()






